In [1]:
from bokeh.plotting import *
import pandas as pd
import numpy as np
from bokeh.layouts import row
df = pd.read_csv('etch_roof_d3s.csv') # The name of the file can be changed depending on the user's needs
arr = df.as_matrix(columns=df.columns[6:]) # Contains all counts per channel
conv = df.as_matrix(columns=df.columns[5:6]).ravel() # Contains list of calibration constants
In [2]:
# arr contains all the counts per channel, but it doesn't include which channel the counts belong to
# In order to account for that, we must create a new array called channels, which is quite simply a list of all the channel numbers
i = 0
channels = []
while i <=1023:
channels.append(i)
i= i+1
######################################################################################################################
# This code aims to sort counts into bins based on the energy of their channel, and in order to do this, we must create an
# array containing all the bin edges
m = 200
bins = []
while m<=2600:
bins.append(m)
m =m+5
######################################################################################################################
buckets = [0] * 1300 # This is the empty array that all our counts will eventually be sorted into
In [3]:
# This cell multiplies the list of channels by each element in the list of calibration factors and sums the results into the previously created energy bins
k = 0
while k < len(conv):
if k%100 == 0:
print k
kev = [i * conv[k] for i in channels]
indexes = np.digitize(kev, bins).flatten()
i = 0
while i < len(indexes):
buckets[indexes[i]-1] = buckets[indexes[i]-1] + arr[k][i]
i = i+1
k = k+1
# Bismuth 214 has a peak at around 620 kev
bismuthA = buckets[70:100]
# -----------------------------------------------------------------------------------------------------------------------
# Potassium 40 has a peak at around 1400 kev
potassiumA = buckets[230:270]
# -----------------------------------------------------------------------------------------------------------------------
# Bismuth 214 has a peak at around 2400 kev
thoriumA = buckets[423:450]
In [4]:
p = figure(y_axis_type="log",title="Entire Spectrum")
p.circle(bins,buckets)
p1 = figure(y_axis_type="log",title="Bismuth Peak")
p1.circle(bins[70:97],bismuthA,color = 'blue')
p2 = figure(y_axis_type="log",title="Potassium Peak")
p2.circle(bins[230:270],potassiumA,color = 'red')
# -----------------------------------------------------------------------------------------------------------------------
# Used to calibrate the thorium peak
# In this code, it is not important
#point = [(bins[423],thoriumA[0]),(bins[450],thoriumA[len(thoriumA)-1])]
#m = ((math.log(18475)-math.log(19395))/(2450-2315))
#c = 10.7061162898
#A = math.exp(c)
#xarino = np.arange(2315,2450)
#y = (A*np.exp((m*xarino)))
#print np.trapz(y,x=xarino)
# -----------------------------------------------------------------------------------------------------------------------
p3 = figure(y_axis_type="log",title="Thallium Peak")
p3.circle(bins[423:450],thoriumA, color = 'green')
show(row(p,p1,p2,p3))